
clear; clc; close all;

%% ========== 1. 系统及分类参数设置 ==========

% (1) 系统参数 (请根据论文实际情况修改)
u = 0.1;      % 控制参数
k = 1.8;      % 耦合强度
a = -1;       % memristor 参数
b = 18;     % memristor 参数
N_total   = 1500;  % 每个初始条件迭代总步数
N_discard = 300;   % 舍弃前 300 步

% (2) 初始条件网格 (请根据论文实际情况修改)
Nx = 200; 
Ny = 200;
x0_min = 0;   x0_max = 0.5;
q0_min = -0.35; q0_max = 0.35;
x0_vals = linspace(x0_min, x0_max, Nx);
q0_vals = linspace(q0_min, q0_max, Ny);

% (3) 分类判定阈值 (可根据系统特性调试)
th_unbounded = 1e5;    % 判断无界的阈值
th_stable    = 1e-2;   % 判断稳定点的 std(x) 阈值
th_period    = 1e-3;   % 周期检测阈值
max_period   = 20;     % 检测的最大周期阶数

%% ========== 2. 扫描 (x0,q0) 并分类 ==========

% 创建存储分类结果的矩阵, 行对应 q0, 列对应 x0
basinData = zeros(Ny, Nx);

for i = 1:Nx
    for j = 1:Ny
        x0_val = x0_vals(i);
        q0_val = q0_vals(j);
        
        % 迭代系统
        [xArr, qArr] = iterate_nhmm(u, k, a, b, x0_val, q0_val, N_total);
        
        % 取后段轨迹
        x_final = xArr(N_discard+1 : end);
        q_final = qArr(N_discard+1 : end);
        
        % 分类
        basinData(j,i) = classifyAttractor(x_final, q_final, ...
                                           th_unbounded, th_stable, ...
                                           th_period, max_period);
    end
    
    % 显示进度
    fprintf('Progress: %.1f%%\n', i/Nx*100);
end

%% ========== 3. 绘制吸引域图 ==========

figure('Position',[100,100,700,600]);
imagesc(x0_vals, q0_vals, basinData);
set(gca,'YDir','normal');  % 让 q0 从下到上
xlabel('x_0'); ylabel('q_0');
title('吸引域分布');

% 定义 4 行 colormap, 顺序与分类值 1~4 对应
% 1=稳定点(粉色), 2=周期(午夜紫), 3=混沌(淡青蓝), 4=无界(棕色)
cmap = [
    1.0  0.8  0.8;  % 稳定点 (Pastel Pink)
    0.2  0.0  0.4;  % 周期 (Midnight Purple)
    0.6  0.9  1.0;  % 混沌 (Light Cyan-Blue)
    0.6  0.3  0.1;  % 无界 (Brown)
];
colormap(cmap);

% 确保 caxis 覆盖 [1,4], 否则最高类可能不显示
caxis([1,4]);

% colorbar
cb = colorbar('Ticks',[1,2,3,4], ...
    'TickLabels',{'稳定点','周期','混沌','无界'});

%% ========== 子函数 2: classifyAttractor ==========
function cls = classifyAttractor(x_final, q_final, th_unbd, th_stable, th_period, max_p)
% 根据最终轨迹 x_final,q_final 分类:
%   1 = 稳定点
%   2 = 周期
%   3 = 混沌
%   4 = 无界
%
% 输入:
%   th_unbd   - 无界判断阈值
%   th_stable - 稳定点判定 std(x) 阈值
%   th_period - 周期检测阈值
%   max_p     - 检测的最大周期阶

% 1) 无界检测
if any(abs(x_final) > th_unbd) || any(abs(q_final) > th_unbd)
    cls = 4; % 无界
    return;
end

% 2) 稳定点: 若 std(x_final) 很小
if std(x_final) < th_stable
    cls = 1; % 稳定点
    return;
end

% 3) 周期检测: 在 2..max_p 范围内尝试
%   如果末段序列与其 p 步延迟之差的最大值 < th_period, 判为周期
is_periodic = false;
len = length(x_final);
for p = 2:max_p
    if len > p
        seg1 = x_final(end-2*p+1 : end-p);
        seg2 = x_final(end-p+1 : end);
        if max(abs(seg1 - seg2)) < th_period
            is_periodic = true;
            break;
        end
    end
end

if is_periodic
    cls = 2; % 周期
else
    cls = 3; % 混沌
end
end